W niniejszej analizie wykorzystano zbiór danych pbc z pakietu Survival. Badany jest wpływ poszczególnych zmiennych, takich jak m.in. wiek, stadium choroby czy związane z badanym schorzeniem symptomy na przeżycie osób cierpiących na pierwotne stwardniające zapalenie dróg żółciowych. Pierwotne stwardniające zapalenie dróg żółciowych (ang. primary sclerosing cholangitis) jest chorobą autoimmunologiczną prowadzącą do wyniszczenia dróg żółciowych w wątrobie. Rozwój choroby jest powolny, lecz nieunikniony, prowadzący do marskości, a nawet wyniszczenia wątroby. Dane pochodzą z badania w Mayo Clinic przeprowadzonego w latach 1974-1984. Jednostką badania jest zbiorowość 424 pacjentów chorujących na PBC skierowanych na badanie pomiędzy 1974 a 1984 rokiem. W analizie przeżycia pod uwagę wzięto tylko pierwsze 312 przypadków, dla których dane są najbardziej kompletne. Zbiór pbc zawiera takie zmienne jak:
Celem badania jest wyodrębnienie czynników wpływających na przeżywalność osób cierpiących na pierwotne stwardniające zapalenie dróg żółciowych.
Biblioteki wykorzystane do realizacji projektu:
library(survival)
library(survminer)
library(survMisc)
library(dplyr)
library(gtools)
library(ggplot2)
library(corrplot)
library(gridExtra)
library(extrafont)
library(VIM)
library(extrafont)
library(mfp)
library(splines)
library(CoxR2)
library(ipred)
library(Rcpp)
library(RColorBrewer)
library(knitr)
library(kableExtra)
library(MASS)
data(pbc)
dane <- pbc
dane$cens <- ifelse(dane$status=="0" | dane$status=="1",0,1)
dane <- dane[1:312,]
Poniżej zaprezentowano 10 pierwszych i 10 ostatnich rekordów ze zbioru danych pbc. W zbiorze występują zmienne zarówno jakościowe, jak i ilościowe. Większość zmiennych jakościowych (wyjątek: sex) zakodowana została za pomocą cyfr “0” i “1” i na tym etapie pracy występują jako zmienne typu numerycznego. Do dalszej części projektu konieczne będzie ich przekształcenie.
head(dane,10)
## id time status trt age sex ascites hepato spiders edema bili chol
## 1 1 400 2 1 58.76523 f 1 1 1 1.0 14.5 261
## 2 2 4500 0 1 56.44627 f 0 1 1 0.0 1.1 302
## 3 3 1012 2 1 70.07255 m 0 0 0 0.5 1.4 176
## 4 4 1925 2 1 54.74059 f 0 1 1 0.5 1.8 244
## 5 5 1504 1 2 38.10541 f 0 1 1 0.0 3.4 279
## 6 6 2503 2 2 66.25873 f 0 1 0 0.0 0.8 248
## 7 7 1832 0 2 55.53457 f 0 1 0 0.0 1.0 322
## 8 8 2466 2 2 53.05681 f 0 0 0 0.0 0.3 280
## 9 9 2400 2 1 42.50787 f 0 0 1 0.0 3.2 562
## 10 10 51 2 2 70.55989 f 1 0 1 1.0 12.6 200
## albumin copper alk.phos ast trig platelet protime stage cens
## 1 2.60 156 1718.0 137.95 172 190 12.2 4 1
## 2 4.14 54 7394.8 113.52 88 221 10.6 3 0
## 3 3.48 210 516.0 96.10 55 151 12.0 4 1
## 4 2.54 64 6121.8 60.63 92 183 10.3 4 1
## 5 3.53 143 671.0 113.15 72 136 10.9 3 0
## 6 3.98 50 944.0 93.00 63 NA 11.0 3 1
## 7 4.09 52 824.0 60.45 213 204 9.7 3 0
## 8 4.00 52 4651.2 28.38 189 373 11.0 3 1
## 9 3.08 79 2276.0 144.15 88 251 11.0 2 1
## 10 2.74 140 918.0 147.25 143 302 11.5 4 1
tail(dane,10)
## id time status trt age sex ascites hepato spiders edema bili chol
## 303 303 1250 0 2 60.65982 f 0 1 1 0 1.0 372
## 304 304 1230 0 1 35.53457 f 0 0 0 0 0.5 219
## 305 305 1216 0 2 43.06639 f 0 1 1 0 2.9 426
## 306 306 1216 0 2 56.39151 f 0 1 0 0 0.6 239
## 307 307 1149 0 2 30.57358 f 0 0 0 0 0.8 273
## 308 308 1153 0 1 61.18275 f 0 1 0 0 0.4 246
## 309 309 994 0 2 58.29979 f 0 0 0 0 0.4 260
## 310 310 939 0 1 62.33265 f 0 0 0 0 1.7 434
## 311 311 839 0 1 37.99863 f 0 0 0 0 2.0 247
## 312 312 788 0 2 33.15264 f 0 0 1 0 6.4 576
## albumin copper alk.phos ast trig platelet protime stage cens
## 303 3.25 108 1190 140 55 248 10.6 4 0
## 304 3.93 22 663 45 75 246 10.8 3 0
## 305 3.61 73 5184 288 144 275 10.6 3 0
## 306 3.45 31 1072 55 64 227 10.7 2 0
## 307 3.56 52 1282 130 59 344 10.5 2 0
## 308 3.58 24 797 91 113 288 10.4 2 0
## 309 2.75 41 1166 70 82 231 10.8 2 0
## 310 3.35 39 1713 171 100 234 10.2 2 0
## 311 3.16 69 1050 117 88 335 10.5 2 0
## 312 3.79 186 2115 136 149 200 10.8 2 0
Sprawdzone zostały także braki danych, których najwięcej jest w przypadku zmiennej dotyczącej poziomu cholesterolu (28 braków) oraz ilości trójglicerydów w organizmie (30 braków). Nieco mniejszą liczbą braków danych obarczone były takie zmienne jak ilość miedzi w moczu (2 braki) oraz ilość płytek krwi (4 braki).
NA_count <- colSums(is.na(dane))
NA_count
## id time status trt age sex ascites hepato
## 0 0 0 0 0 0 0 0
## spiders edema bili chol albumin copper alk.phos ast
## 0 0 0 28 0 2 0 0
## trig platelet protime stage cens
## 30 4 0 0 0
barplot(NA_count)
W celu pozbycia się braków danych wykorzystano funkcję kNN z pakietu VIM – imputację danych wykonano za pomocą metody k-najbliższych sąsiadów (parametr k ustawiono jako 5). Co więcej, zmienne wcześniej zakodowane jako zmienne ilościowe, przekształcono w zmienne jakościowe (status, trt, ascites, spiders, hepato, stage, edema, cens).
daneKNN <- kNN(dane, variable = c("chol", "copper", "trig", "platelet"), k = 5)
daneKNN <- daneKNN[,1:21,drop = F]
daneKNN$age <- round(daneKNN$age,0)
levels(daneKNN$sex) <- c(1,0)
daneKNN$status <- as.factor(daneKNN$status)
daneKNN$trt <- as.factor(daneKNN$trt)
daneKNN$ascites <- as.factor(daneKNN$ascites)
daneKNN$spiders <- as.factor(daneKNN$spiders)
daneKNN$hepato <- as.factor(daneKNN$hepato)
daneKNN$stage <- as.factor(daneKNN$stage)
daneKNN$edema <- as.factor(daneKNN$edema)
daneKNN$cens <- as.factor(daneKNN$cens)
str(daneKNN)
## 'data.frame': 312 obs. of 21 variables:
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ time : int 400 4500 1012 1925 1504 2503 1832 2466 2400 51 ...
## $ status : Factor w/ 3 levels "0","1","2": 3 1 3 3 2 3 1 3 3 3 ...
## $ trt : Factor w/ 2 levels "1","2": 1 1 1 1 2 2 2 2 1 2 ...
## $ age : num 59 56 70 55 38 66 56 53 43 71 ...
## $ sex : Factor w/ 2 levels "1","0": 2 2 1 2 2 2 2 2 2 2 ...
## $ ascites : Factor w/ 2 levels "0","1": 2 1 1 1 1 1 1 1 1 2 ...
## $ hepato : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 2 1 1 1 ...
## $ spiders : Factor w/ 2 levels "0","1": 2 2 1 2 2 1 1 1 2 2 ...
## $ edema : Factor w/ 3 levels "0","0.5","1": 3 1 2 2 1 1 1 1 1 3 ...
## $ bili : num 14.5 1.1 1.4 1.8 3.4 0.8 1 0.3 3.2 12.6 ...
## $ chol : int 261 302 176 244 279 248 322 280 562 200 ...
## $ albumin : num 2.6 4.14 3.48 2.54 3.53 3.98 4.09 4 3.08 2.74 ...
## $ copper : int 156 54 210 64 143 50 52 52 79 140 ...
## $ alk.phos: num 1718 7395 516 6122 671 ...
## $ ast : num 137.9 113.5 96.1 60.6 113.2 ...
## $ trig : int 172 88 55 92 72 63 213 189 88 143 ...
## $ platelet: int 190 221 151 183 136 233 204 373 251 302 ...
## $ protime : num 12.2 10.6 12 10.3 10.9 11 9.7 11 11 11.5 ...
## $ stage : Factor w/ 4 levels "1","2","3","4": 4 3 4 4 3 3 3 3 2 4 ...
## $ cens : Factor w/ 2 levels "0","1": 2 1 2 2 1 2 1 2 2 2 ...
summary(daneKNN)
## id time status trt age sex
## Min. : 1.00 Min. : 41 0:168 1:158 Min. :26.00 1: 36
## 1st Qu.: 78.75 1st Qu.:1191 1: 19 2:154 1st Qu.:42.00 0:276
## Median :156.50 Median :1840 2:125 Median :50.00
## Mean :156.50 Mean :2006 Mean :50.03
## 3rd Qu.:234.25 3rd Qu.:2697 3rd Qu.:57.00
## Max. :312.00 Max. :4556 Max. :78.00
## ascites hepato spiders edema bili chol
## 0:288 0:152 0:222 0 :263 Min. : 0.300 Min. : 120.0
## 1: 24 1:160 1: 90 0.5: 29 1st Qu.: 0.800 1st Qu.: 252.0
## 1 : 20 Median : 1.350 Median : 309.0
## Mean : 3.256 Mean : 364.8
## 3rd Qu.: 3.425 3rd Qu.: 394.2
## Max. :28.000 Max. :1775.0
## albumin copper alk.phos ast
## Min. :1.96 Min. : 4.00 Min. : 289.0 Min. : 26.35
## 1st Qu.:3.31 1st Qu.: 41.75 1st Qu.: 871.5 1st Qu.: 80.60
## Median :3.55 Median : 73.00 Median : 1259.0 Median :114.70
## Mean :3.52 Mean : 97.73 Mean : 1982.7 Mean :122.56
## 3rd Qu.:3.80 3rd Qu.:123.25 3rd Qu.: 1980.0 3rd Qu.:151.90
## Max. :4.64 Max. :588.00 Max. :13862.4 Max. :457.25
## trig platelet protime stage cens
## Min. : 33.00 Min. : 62.0 Min. : 9.00 1: 16 0:187
## 1st Qu.: 84.75 1st Qu.:200.0 1st Qu.:10.00 2: 67 1:125
## Median :107.50 Median :256.0 Median :10.60 3:120
## Mean :122.34 Mean :261.9 Mean :10.73 4:109
## 3rd Qu.:146.00 3rd Qu.:322.5 3rd Qu.:11.10
## Max. :598.00 Max. :563.0 Max. :17.10
a <- ggplot(data = daneKNN, aes(x = age)) + geom_histogram(alpha = 0.7, binwidth = 5, fill="darkorange") + theme_minimal() +
ggtitle('Age') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
b <- ggplot(data = daneKNN, aes(x = bili)) + geom_histogram(alpha = 0.7, binwidth = 2, fill="mediumpurple3") + theme_minimal() +
ggtitle('Bilirunbin') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
c <- ggplot(data = daneKNN, aes(x = chol)) + geom_histogram(alpha = 0.7, binwidth = 50, fill="seagreen4") + theme_minimal() +
ggtitle('Cholesterol') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
d <- ggplot(data = daneKNN, aes(x = albumin)) + geom_histogram(alpha = 0.7, binwidth = 0.2, fill="dodgerblue3") + theme_minimal() +
ggtitle('Albumin') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
a1 <- ggplot(data = daneKNN, aes(x = copper)) + geom_histogram(alpha = 0.7, binwidth = 30, fill="indianred4") + theme_minimal() +
ggtitle('Copper') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "N")
b1 <- ggplot(data = daneKNN, aes(x = alk.phos)) + geom_histogram(alpha = 0.7, binwidth = 400, fill="royalblue3") + theme_minimal() +
ggtitle('Alkaline phos.') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
c1 <- ggplot(data = daneKNN, aes(x = ast)) + geom_histogram(alpha = 0.7, binwidth = 20, fill="goldenrod2") + theme_minimal() +
ggtitle('Aspartate amino.') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
d1 <- ggplot(data = daneKNN, aes(x = trig)) + geom_histogram(alpha = 0.7, binwidth = 20, fill="cadetblue3") + theme_minimal() +
ggtitle('Triglycerides') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
a2 <- ggplot(data = daneKNN, aes(x = platelet)) + geom_histogram(alpha = 0.7, binwidth = 25, fill="orangered3") + theme_minimal() +
ggtitle('Platelet') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
b2 <- ggplot(data = daneKNN, aes(x = protime)) + geom_histogram(alpha = 0.7, binwidth = 0.5, fill="darkolivegreen4") + theme_minimal() +
ggtitle('Protime') +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
c2 <- ggplot(data = daneKNN, aes(x = time, fill = cens)) + geom_histogram(alpha = 0.7, binwidth = 200) + facet_grid(~cens) + theme_minimal() +
ggtitle('Time for cens & event') +
scale_fill_manual(values = c("coral3","slateblue3")) +
theme(axis.title.y = element_text(color="Grey23", size=11),
axis.text.y = element_text(size=9),
axis.title.x = element_blank(),
plot.title = element_text(color="gray26", size=20, family="serif")) +
labs(y = "N")
W przypadku zmiennej dotyczącej wieku, najniższą wartością jest 26, najwyższą zaś 78. Średni wiek badanych osób to 50 lat i wokół tej wartości skupia się większość obserwacji. Poziom bilirubiny w normie waha się od 0,2 mg/dl do 1,1 mg/dl i w przypadku połowy obserwacji jej wartości nie przekraczają 1,35 mg/dl. Natomiast obserwacja, u której występuje najwyższa wartość bilirubiny to aż 28 mg/dl. Podobnie jest w przypadku cholesterolu, którego norma powinna wynosić mniej niż 200 mg/dl. Pomimo że rozkład zmiennej charakteryzuje się silną asymetrią prawostronną, to u większości jednostek cholesterol przekracza zakładaną normę, bowiem wartość pierwszego kwartyla to już 252 mg/dl. Wartość maksymalna przekracza normę niemal dziewięciokrotnie (1775 mg/dl). Zmienna albumin posiada rozkład lewostronny, choć w jego przypadku asymetria nie jest duża.
grid.arrange(a,b,c,d)
Zmienne widoczne na poniższych wykresach charakterysują sie asymetrą prawostronną. Najmniejsza asymetria prawostronna widoczna jest dla zmiennej aspartate amino..
grid.arrange(a1,b1,c1,d1)
Z histogramu dotyczącego zmiennej platelet wynika że w przypadku większości pacjentów ilość płytek krwi (trombocytów) była w normie (150 – 400 tys./mikrolitr). Podobnie jest w przypadku zmiennej protime, której 50% wartości mieściło się w przedziale [10, 11,1], a normą czasu krzepnięcia krwi w teście zastosowanym w badaniu Mayo jest przedział [9,5, 11,8] sekundy.
grid.arrange(a2,b2)
Rozkłady zmiennych dotyczących zarówno cenzurowania jak i zdarzenia względem czasu trwania nie różnią się znacząco.
c2
e <- daneKNN %>%
count(sex) %>%
ggplot() + geom_col(aes(x = sex, y = n, fill = sex), alpha = 0.7)+ geom_label(aes(x = sex, y = n, label = n)) + scale_fill_manual(values = c("navajowhite3","plum4")) +
theme_minimal() +
ggtitle('Number of patients by gender') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
f <- daneKNN %>%
count(trt) %>%
ggplot() + geom_col(aes(x = trt, y = n, fill = trt), alpha = 0.7)+ geom_label(aes(x = trt, y = n, label = n)) + scale_fill_manual(values = c("lightseagreen","cornsilk3")) +
theme_minimal() +
ggtitle('Number of patients by treatment') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
g <- daneKNN %>%
count(ascites) %>%
ggplot() + geom_col(aes(x = ascites, y = n, fill = ascites), alpha = 0.7)+ geom_label(aes(x = ascites, y = n, label = n)) + scale_fill_manual(values = c("darkseagreen4","coral3")) +
theme_minimal() +
ggtitle('Number of patients by ascites presence') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
h <- daneKNN %>%
count(hepato) %>%
ggplot() + geom_col(aes(x = hepato, y = n, fill = hepato), alpha = 0.7)+ geom_label(aes(x = hepato, y = n, label = n)) + scale_fill_manual(values = c("palegreen4", "indianred3")) +
theme_minimal() +
ggtitle('Number of patients by hepato presence') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
e1 <- daneKNN %>%
count(edema) %>%
ggplot() + geom_col(aes(x = edema, y = n, fill = edema), alpha = 0.7)+ geom_label(aes(x = edema, y = n, label = n)) + scale_fill_manual(values = c("springgreen3","tan3","tomato3")) +
theme_minimal() +
ggtitle('Number of patients by edema presence') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
f1 <- daneKNN %>%
count(stage) %>%
ggplot() + geom_col(aes(x = stage, y = n, fill = stage), alpha = 0.7)+ geom_label(aes(x = stage, y = n, label = n)) + scale_fill_manual(values = c("palegreen3", "skyblue3","salmon3","brown3")) +
theme_minimal() +
ggtitle('Number of patients by stage of disease') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
g1 <- daneKNN %>%
count(spiders) %>%
ggplot() + geom_col(aes(x = spiders, y = n, fill = spiders), alpha = 0.7)+ geom_label(aes(x = spiders, y = n, label = n)) + scale_fill_manual(values = c("forestgreen", "chocolate3")) +
theme_minimal() +
ggtitle('Number of patients spiders presence') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
h1 <- daneKNN %>%
count(cens) %>%
ggplot() + geom_col(aes(x = cens, y = n, fill = cens), alpha = 0.7)+ geom_label(aes(x = cens, y = n, label = n)) + scale_fill_manual(values = c("coral3","slateblue3")) +
theme_minimal() +
ggtitle('Number of patients by cens or event') +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=19, family="serif"))
Jeśli chodzi o rozkłady zmiennych jakościowych, obserwuje się znacznie większą liczbę pacjentów płci żeńskiej (276 jednostek) w porównaniu do liczby pacjentów płci męskiej (zaledwie 36 jednostek). Niemal połowa badanych jednostek otrzymała placebo (154 osoby), 158 osób zostało poddanych działaniu leku D-penicillmain.
grid.arrange(e,f)
Spośród 312 badanych jednostek tylko 24 osoby cierpiały na wodobrzusze, natomiast jeżeli chodzi o obecność powiększonej wątroby, choroba ta została wykryta już u 160 pacjentów.
grid.arrange(g,h)
U większości badanych (263 jednostki) obrzęk spowodowany nadmiarem płynu uwięzionego w tkankach organizmu nie występuje. W przypadku 29 osób schorzenie zostało wyleczone lub było nieleczone, a u 20 osób edema występuje pomimo zastosowania terapii diuretykami. Najwięcej badanych znajduje się w późniejszych stadiach choroby: III (120 jednostek) oraz IV (109 jednostek). Najmniej jednostek znajduje się w I stadium choroby i jest to 16 osób spośród wszystkich badanych.
grid.arrange(e1,f1)
Stosunkowo więcej jednostek niż w przypadku wyżej wymienionych schorzeń oraz objawów współwystępujących charakteryzuje się występowaniem zmian naczyniowych, bowiem to już 90 badanych. 125 spośród badanych jednostek to obserwacje cenzurowane.
grid.arrange(g1,h1)
z <- ggplot(daneKNN, aes(x=stage, y = age, fill=stage)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen3", "skyblue3","salmon3","brown3")) + theme_minimal() +
ggtitle("Age by stage") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Age")
x <- ggplot(daneKNN, aes(x=stage, y = bili, fill=stage)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen3", "skyblue3","salmon3","brown3")) + theme_minimal() +
ggtitle("Bilirunbin by stage") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Bilirunbin")
z1 <- ggplot(daneKNN, aes(x=stage, y = time, fill=stage)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen3", "skyblue3","salmon3","brown3")) + theme_minimal() +
ggtitle("Time by stage") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Time") + facet_grid(~cens)
x1 <- ggplot(daneKNN, aes(x=ascites, y = bili, fill=ascites)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("darkseagreen4","coral3")) + theme_minimal() +
ggtitle("Bilirunbin by ascites presence") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Bilirunbin")
x2 <- ggplot(daneKNN, aes(x=hepato, y = bili, fill=hepato)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen4", "indianred3")) + theme_minimal() +
ggtitle("Bilirunbin by hepato presence") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Bilirunbin")
z2 <- ggplot(daneKNN, aes(x=hepato, y = copper, fill=hepato)) + geom_boxplot(size=1.2, alpha=0.5) +
scale_fill_manual(values = c("palegreen4", "indianred3")) + theme_minimal() +
ggtitle("Copper by hepato presence") +
theme(axis.title.y = element_text(color="Grey23", size=12),
axis.title.x=element_blank(),
legend.title = element_text(size=12),
legend.text = element_text(size=8),
legend.position = "right",
legend.justification = c(0.94,0.94),
legend.background = element_rect(fill="grey88",
size=0.5, linetype="solid",
colour ="gray26"),
plot.title = element_text(color="gray26", size=18, family="serif")) +
labs(y = "Copper")
Najmłodsze z badanych osób znajdują się raczej w I stadium choroby, najstarsze zaś w ostatnim IV. Najwyższymi wartościami poziomu barwnika żółciowego charakteryzują się osoby z ostatniego stadia choroby.
grid.arrange(z,x)
W przypadku liczby dni pomiędzy rejestracją a śmiercią, transplantacją lub innym powodem przerwania badania – im wyższe stadium choroby, tym krótszy czas obserwacji.
Osoby z wodobrzuszem charakteryzują się dużo wyższymi wartościami stężenia barwnika żółciowego we krwi.
grid.arrange(z1,x1)
Podobna zależność występuje w przypadku związku ilości miedzi w moczu, a występowaniem powiększenia wątroby. Analogicznie jest również w przypadku zależności pomiędzy stężeniem barwnika żółciowego, a występowaniem powiększenia wątroby.
grid.arrange(z2,x2)
Jeśli chodzi o różnice w wieku pomiędzy grupami kobiet i mężczyzn, nie są one wysokie, w grupie mężczyzn znajduje się jednostka najstarsza, natomiast w grupie kobiet jednostka najmłodsza. Średnia różni się pomiędzy grupami o 7 lat.
tab1 <- daneKNN %>%
group_by(sex) %>%
summarise(N = n(), Mean = round(mean(age),2), Min = min(age), Max = max(age), Sd = round(sd(age),2))
tab1 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)
| sex | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 36 | 56.19 | 33 | 78 | 11.54 |
| 0 | 276 | 49.22 | 26 | 77 | 10.20 |
Stężenie barwnika żółciowego jest średnio wyższe u osób, u których występuje obrzęk spowodowany nadmiarem płynu uwięzionego w tkankach organizmu. Co więcej, odchylenie standardowe wskazuje na to, że wartości stężenia barwnika żółciowego u osób, u których edema nie występuje, nie różnią się tak bardzo wewnątrz tej grupy, jak w przypadku grup, u których schorzenie zostało wyleczone lub było nieleczone czy też u których edema występuje pomimo zastosowania terapii diuretykami.Stężenie barwnika żółciowego jest średnio wyższe u osób, u których występuje obrzęk spowodowany nadmiarem płynu uwięzionego w tkankach organizmu. Co więcej, odchylenie standardowe wskazuje na to, że wartości stężenia barwnika żółciowego u osób, u których edema nie występuje, nie różnią się tak bardzo wewnątrz tej grupy, jak w przypadku grup, u których schorzenie zostało wyleczone lub było nieleczone czy też u których edema występuje pomimo zastosowania terapii diuretykami.
tab2 <- daneKNN %>%
group_by(edema) %>%
summarise(N = n(), Mean = round(mean(bili),2), Min = min(bili), Max = max(bili), Sd = round(sd(bili),2))
tab2 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)
| edema | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 0 | 263 | 2.52 | 0.3 | 25.5 | 3.35 |
| 0.5 | 29 | 5.76 | 0.6 | 28.0 | 7.24 |
| 1 | 20 | 9.26 | 0.8 | 22.5 | 7.00 |
Średnie stężenie barwnika żółciowego jest największe u osób, które posiadają najwyższe stadium rozwoju choroby.
tab3 <- daneKNN %>%
group_by(stage) %>%
summarise(N = n(), Mean = round(mean(bili),2), Min = min(bili), Max = max(bili), Sd = round(sd(bili),2))
tab3 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)
| stage | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 16 | 1.07 | 0.5 | 6.0 | 1.33 |
| 2 | 67 | 2.13 | 0.3 | 25.5 | 3.82 |
| 3 | 120 | 2.90 | 0.3 | 28.0 | 4.34 |
| 4 | 109 | 4.66 | 0.6 | 24.5 | 5.05 |
Kobiety charakteryzują się stosunkowo nieco wyższym poziomem cholesterolu niż mężczyźni, choć może to wynikać z dysproporcji w liczebności obu grup, bowiem również wartości skrajne są dużo wyższe.
tab4 <- daneKNN %>%
group_by(sex) %>%
summarise(N = n(), Mean = round(mean(chol),2), Min = min(chol), Max = max(chol), Sd = round(sd(chol),2))
tab4 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)
| sex | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 36 | 357.61 | 151 | 1000 | 178.80 |
| 0 | 276 | 365.79 | 120 | 1775 | 228.08 |
Jeśli chodzi o ilość trójglicerydów względem płci, średnia jest wyższa w przypadku mężczyzn, mimo że np. maksymalna wartość dla mężczyzn wynosi 242 mg/dl, a u kobiet jest dwukrotnie wyższa – 598 mg/dl.
tab5 <- daneKNN %>%
group_by(sex) %>%
summarise(N = n(), Mean = round(mean(trig),2), Min = min(trig), Max = max(trig), Sd = round(sd(trig),2))
tab5 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)
| sex | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 36 | 132.56 | 49 | 242 | 51.68 |
| 0 | 276 | 121.00 | 33 | 598 | 64.12 |
Średnie stężenie trójglicerydów we krwi rośnie również wraz ze stadium rozwoju choroby, choć średnia w III stadium jest nieco wyższa niż w stadium IV.
tab6 <- daneKNN %>%
group_by(stage) %>%
summarise(N = n(), Mean = round(mean(trig),2), Min = min(trig), Max = max(trig), Sd = round(sd(trig),2))
tab6 %>%
kbl() %>%
kable_styling(full_width = T, html_font = "Cambria", fixed_thead = T, font_size = 16)
| stage | N | Mean | Min | Max | Sd |
|---|---|---|---|---|---|
| 1 | 16 | 90.75 | 55 | 188 | 32.22 |
| 2 | 67 | 112.63 | 44 | 319 | 53.19 |
| 3 | 120 | 128.68 | 46 | 382 | 56.91 |
| 4 | 109 | 125.96 | 33 | 598 | 75.28 |
Najsilniejsze korelacje występują w przypadku takich par zmiennych jak: protime, copper (0,22); ast, copper (0,3); trig, copper (0,29); ast, chol (0,35); trig, chol (0,28) czy też protime, platelet (-0,22); protime, albumin (-0,23); copper, albumin (-0,27); ast, albumin (-0,22).
doKor <- daneKNN[,c(5, 12,13,14,15,16,17,18,19)]
corMat <- cor(doKor)
#corrplot(corMat,method="number", tl.col = 'black')
corrplot(corMat, method = 'color',addCoef.col = 'black', order = 'AOE', tl.col = 'black', col=brewer.pal(n=10, name="PuOr"))
Na wykresach zostały pokazane funkcje przeżycia dla wszystkich jednostek z wykorzystaniem estymatora Kaplana - Meiera oraz Nelsona - Aelena. Dodatkowo, do oszacowania funkcji przeżycia za pomocą estymatora KM zostały wykorzystane dwa rodzaje przedziałów ufności - liniowy i logarytmiczny. Na wykresach brak jest istotnych różnic pomiędzy porównywanymi estymatorami. Potwierdza to również średnia oraz mediana czasu trwania widoczna na podsumowaniu pod wykresami. W momencie w którym wykresy się urywają (po 4000 dni) pozostaje już tylko około 25% jednostek.
daneKNN$cens <- as.character(daneKNN$cens)
daneKNN$cens <- as.numeric(daneKNN$cens)
pbcKM <- survfit(Surv(time, cens) ~ 1, conf.type="plain", data=daneKNN) #KM
pbcKM1 <- survfit(Surv(time, cens) ~ 1, conf.type="log-log", data=daneKNN)
pbcNELS=survfit(Surv(time,cens) ~ 1, conf.type ="plain", type="fleming-harrington",data=daneKNN)
splots <- list()
splots[[1]] <- ggsurvplot(
fit = pbcKM,
xlab = "Days",
ylab = "Overall survival probability with Kaplan - Meier estimator and plain confidence interval")
splots[[2]] <- ggsurvplot(
fit = pbcKM1,
xlab = "Days",
ylab = "Overall survival probability with with Kaplan - Meier estimator and log confidence interval",
palette = "blue")
splots[[3]] <- ggsurvplot(
fit = pbcNELS,
xlab = "Days",
ylab = "Overall survival probability with Nelson - Aelen estimator",
palette = "limegreen")
arrange_ggsurvplots(splots, print = TRUE,
ncol = 3, nrow = 1)
summary(pbcKM)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL
## 312.0000 312.0000 312.0000 125.0000 2973.6112 101.8206 3395.0000 3086.0000
## 0.95UCL
## 3839.0000
summary(pbcKM1)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL
## 312.0000 312.0000 312.0000 125.0000 2973.6112 101.8206 3395.0000 3086.0000
## 0.95UCL
## 3839.0000
summary(pbcNELS)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL
## 312.0000 312.0000 312.0000 125.0000 2978.9327 102.2104 3395.0000 3086.0000
## 0.95UCL
## 3839.0000
daneKNN$Age_cat=quantcut(daneKNN$age,q=4)
daneKNN$bili_cat=quantcut(daneKNN$bili,q=4)
daneKNN$chol_cat=quantcut(daneKNN$chol,q=4)
daneKNN$albumin_cat=quantcut(daneKNN$albumin,q=4)
daneKNN$copper_cat=quantcut(daneKNN$copper,q=4)
daneKNN$alk_cat=quantcut(daneKNN$alk.phos,q=4)
daneKNN$ast_cat=quantcut(daneKNN$ast,q=4)
daneKNN$trig_cat=quantcut(daneKNN$trig,q=4)
daneKNN$plat_cat=quantcut(daneKNN$platelet,q=4)
daneKNN$pro_cat=quantcut(daneKNN$protime,q=4)
Grupy w zmiennych jakościowych zostały wydorębnione na podstawie wariantów zakodowanych w zmienych, natomiast zmienne ilościowe zostały podzielone za pomocą funkcji quantcut na 4 przedziały kwartylowe i każdy przedział stanowił odrębny wariant nowopowstałej zmiennej.
Dla zmiennej sex większe prawdopodobieństwo przeżycia niemal w każdym momencie trwania jest w grupie kobiet. Mediana czasu przeżycia dla mężczyzn wynosi 2386 dni, natomiast dla kobiet 3428 dni. Test Log-Rank wskazał na istotne różnice pomiędzy tymi grupami jednostek, więc rozkłady czasu trwania są różne (p-value = 0.039). Z kolei dla pacjentów podzielonych względem zmiennej trt, test Log-Rank nie wskazał istotnych różnic (p-value = 0.75).
groupSplots1 <- list()
groupSplots1[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ sex, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("navajowhite3","plum4"))
groupSplots1[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ trt, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("lightseagreen","cornsilk3"))
arrange_ggsurvplots(groupSplots1, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmsex <- survfit(Surv(time, cens) ~ sex, data = daneKNN)
summary(kmsex)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
## sex=1 36 36 36 22 2479.267 281.2595 2386 1297 NA
## sex=0 276 276 276 103 3048.694 109.5150 3428 3170 NA
kmtrt <- survfit(Surv(time, cens) ~ trt, data = daneKNN)
summary(kmtrt)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
## trt=1 158 158 158 65 2949.313 141.7174 3282 2583 NA
## trt=2 154 154 154 60 3002.749 145.7727 3428 3090 NA
Pacjenci, u których stwierdzono wodobrzusze (ascites) posiadają dużo niższe prawdopodobieństwo przeżycia w każdym momencie trwania. Istnieją również istotne różnice w grupach wydzielonych na podstawie tej zmiennej (p-value < 0.0001). Podobna sytuacja klaruje się w przypadku zmiennej spiders (p-value < 0.0001). Mediana czasu trwania dla pacjentów, u których wykryto malformacje naczyniową wynosi zaledwie 1847 dni.
groupSplots2 <- list()
groupSplots2[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ ascites, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("darkseagreen4","coral3"))
groupSplots2[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ spiders, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("forestgreen", "chocolate3"))
arrange_ggsurvplots(groupSplots2, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmascites <- survfit(Surv(time, cens) ~ ascites, data = daneKNN)
summary(kmascites)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL
## ascites=0 288 288 288 102 3165.2653 102.3598 3584 3282
## ascites=1 24 24 24 23 856.1944 201.0002 368 223
## 0.95UCL
## ascites=0 NA
## ascites=1 1191
kmspiders <- survfit(Surv(time, cens) ~ spiders, data = daneKNN)
summary(kmspiders)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL
## spiders=0 222 222 222 73 3287.546 113.0711 3839 3282
## spiders=1 90 90 90 52 2153.799 188.2438 1847 1434
## 0.95UCL
## spiders=0 NA
## spiders=1 3244
Pacjenci, u których wykryto powiększoną wątrobę (hepato) i obrzęk spowodowany nadmiarem płynu uwięzionego w tkankach organizmu (edema) posiadają dużo niższe prawdopodobieństwo przeżycia w każdym momencie trwania. Również średnia dni trwania w poszczególnych grupach jest dużo niższa dla osób zmagających się z tymi chorobami. Test Log-Rank wykazał na istotne różnice pomiędzy grupami, a ich rozkładami czasu trwania (p-value < 0.0001).
groupSplots3 <- list()
groupSplots3[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ hepato, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("palegreen4", "indianred3"))
groupSplots3[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ edema, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("springgreen3","tan3","tomato3"))
arrange_ggsurvplots(groupSplots3, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmhepato <- survfit(Surv(time, cens) ~ hepato, data = daneKNN)
summary(kmhepato)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
## hepato=0 152 152 152 37 3606.189 129.0008 NA 3584 NA
## hepato=1 160 160 160 88 2372.598 138.0432 2256 1741 3244
kmedema <- survfit(Surv(time, cens) ~ edema, data = daneKNN)
summary(kmedema)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL
## edema=0 263 263 263 89 3231.814 104.2501 3584 3244
## edema=0.5 29 29 29 17 2269.608 350.2181 1576 1012
## edema=1 20 20 20 19 673.550 207.7994 299 179
## 0.95UCL
## edema=0 NA
## edema=0.5 NA
## edema=1 971
Wraz ze wzrostem stadium choroby spada prawdopodobieństwo przeżycia jednostek. Średnia czasu trwania dla pacjentów zmniejsza się, im stadium rozwoju choroby jest większe. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value < 0.0001).
ggsurvplot(survfit(Surv(time, cens) ~ stage, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_bw(),
palette = c("palegreen3", "skyblue3","salmon3","brown3"))
kmstage <- survfit(Surv(time, cens) ~ stage, data = daneKNN)
summary(kmstage)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
## stage=1 16 16 16 1 4358.727 188.0922 NA NA NA
## stage=2 67 67 67 16 3636.000 179.4216 4079 3839 NA
## stage=3 120 120 120 43 3156.503 153.5608 3428 2847 NA
## stage=4 109 109 109 65 2101.587 173.8259 1682 1165 2796
W przypadku zmiennej trig prawdopodobieństwo przetrwania jest większe w grupach, dla których ilość trójglicerydów w organizmie jest niższa. Analogicznie jest w przypadku wieku - im niższy przedział wiekowy, tym większe jest prawdopodobieństwo przeżycia w danej grupie. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value = 0.0024 dla zmiennej trig oraz p-value < 0.0001 dla zmiennej age).
groupSplots4 <- list()
groupSplots4[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ trig_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots4[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ Age_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots4, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmtrig <- survfit(Surv(time, cens) ~ trig_cat, data = daneKNN)
summary(kmtrig)$table
## records n.max n.start events rmean se(rmean) median
## trig_cat=[33,84.8] 78 78 78 22 3515.323 181.6379 NA
## trig_cat=(84.8,108] 78 78 78 32 2968.504 195.0771 3282
## trig_cat=(108,146] 81 81 81 30 2941.081 222.9211 3170
## trig_cat=(146,598] 75 75 75 41 2468.393 203.3881 2466
## 0.95LCL 0.95UCL
## trig_cat=[33,84.8] 3428 NA
## trig_cat=(84.8,108] 2796 NA
## trig_cat=(108,146] 3090 NA
## trig_cat=(146,598] 1682 4079
kmage <- survfit(Surv(time, cens) ~ Age_cat, data = daneKNN)
summary(kmage)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL
## Age_cat=[26,42] 79 79 79 19 3492.100 199.5027 NA 3428
## Age_cat=(42,50] 85 85 85 28 3335.602 184.0169 4191 3358
## Age_cat=(50,57] 76 76 76 36 2831.014 192.2936 3282 2466
## Age_cat=(57,78] 72 72 72 42 2190.896 206.2281 2090 1356
## 0.95UCL
## Age_cat=[26,42] NA
## Age_cat=(42,50] NA
## Age_cat=(50,57] NA
## Age_cat=(57,78] 3222
W przypadku zmiennej bili prawdopodobieństwo przetrwania jest większe w grupach, dla których stężenie barwnika żółciowego w organizmie jest niższe. Podobnie jest w przypadku cholesterolu - im niższa jest wartość zmiennej chol, tym większe jest prawdopodobieństwo przeżycia w danej grupie, z wyjątkiem grupy, dla której cholesterol jest najniższy - w tym wypadku funkcja trwania pokrywa się z pozostałymi. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value < 0.0001 dla zmiennej bili oraz p-value = 0.00067 dla zmiennej chol).
groupSplots5 <- list()
groupSplots5[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ bili_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots5[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ chol_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots5, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmbili <- survfit(Surv(time, cens) ~ bili_cat, data = daneKNN)
summary(kmbili)$table
## records n.max n.start events rmean se(rmean) median
## bili_cat=[0.3,0.8] 90 90 90 14 3964.305 139.3599 NA
## bili_cat=(0.8,1.35] 66 66 66 14 3722.350 188.9378 NA
## bili_cat=(1.35,3.42] 78 78 78 39 2683.480 180.4503 3222
## bili_cat=(3.42,28] 78 78 78 58 1447.146 144.4648 1170
## 0.95LCL 0.95UCL
## bili_cat=[0.3,0.8] NA NA
## bili_cat=(0.8,1.35] 3090 NA
## bili_cat=(1.35,3.42] 2105 3445
## bili_cat=(3.42,28] 943 1427
kmchol <- survfit(Surv(time, cens) ~ chol_cat, data = daneKNN)
summary(kmchol)$table
## records n.max n.start events rmean se(rmean) median
## chol_cat=[120,252] 79 79 79 28 3074.415 220.9585 NA
## chol_cat=(252,309] 78 78 78 24 3465.298 170.6644 3762
## chol_cat=(309,394] 77 77 77 28 2947.261 220.5058 2847
## chol_cat=(394,1.78e+03] 78 78 78 45 2411.435 180.8935 2105
## 0.95LCL 0.95UCL
## chol_cat=[120,252] 3170 NA
## chol_cat=(252,309] 3395 NA
## chol_cat=(309,394] 2419 NA
## chol_cat=(394,1.78e+03] 1492 3574
W przypadku skategoryzowanej zmiennej albumin, najniższe prawdopodobieństwo przeżycia jest w grupie osób posiadających niski poziom wartości tej zmiennej. Test Log-Rank wykazał istotne różnice pomiędzy grupami w zależności od ich rozkładów czasu trwania (p-value < 0.0001). Podobna sytuacja występuje w przypadku zmiennej copper, jeśli chodzi o różne rozkłady czasu trwania (p-value < 0.0001). Jednakże dla tej zmiennej, niskie wartości oznaczają większe prawdopodobieństwo przeżycia.
groupSplots6 <- list()
groupSplots6[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ albumin_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots6[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ copper_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots6, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmalbumin <- survfit(Surv(time, cens) ~ albumin_cat, data = daneKNN)
summary(kmalbumin)$table
## records n.max n.start events rmean se(rmean) median
## albumin_cat=[1.96,3.31] 81 81 81 55 1543.119 131.0277 1217
## albumin_cat=(3.31,3.55] 76 76 76 30 2985.155 203.5664 3358
## albumin_cat=(3.55,3.8] 81 81 81 18 3691.679 171.3645 NA
## albumin_cat=(3.8,4.64] 74 74 74 22 3495.126 170.9366 3762
## 0.95LCL 0.95UCL
## albumin_cat=[1.96,3.31] 1000 2256
## albumin_cat=(3.31,3.55] 2583 NA
## albumin_cat=(3.55,3.8] 3853 NA
## albumin_cat=(3.8,4.64] 3574 NA
kmcopper <- survfit(Surv(time, cens) ~ copper_cat, data = daneKNN)
summary(kmcopper)$table
## records n.max n.start events rmean se(rmean) median
## copper_cat=[4,41.8] 78 78 78 11 4043.353 140.0845 NA
## copper_cat=(41.8,73] 81 81 81 25 3250.959 191.9466 3358
## copper_cat=(73,123] 75 75 75 32 2875.992 205.1588 3170
## copper_cat=(123,588] 78 78 78 57 1734.077 172.5699 1170
## 0.95LCL 0.95UCL
## copper_cat=[4,41.8] NA NA
## copper_cat=(41.8,73] 3090 NA
## copper_cat=(73,123] 2400 NA
## copper_cat=(123,588] 930 1827
Skategoryzowane zmienne alk i ast wykazały, że jednostki należące do poszczególnych grup wyodrębnionych z uwagi na wartości tych zmiennych posiadają różne rozkładu czasu trwania. P-value w teście Log-Rank wyniosło odpowiednio 0.002 i < 0.0001. Pacjenci, którzy posiadają mniejsze wartości zmiennej alk i ast posiadają większe prawdopodobieństwo przeżycia, niż pozostali.
groupSplots7 <- list()
groupSplots7[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ alk_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots7[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ ast_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots7, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmalk <- survfit(Surv(time, cens) ~ alk_cat, data = daneKNN)
summary(kmalk)$table
## records n.max n.start events rmean se(rmean)
## alk_cat=[289,872] 78 78 78 18 3589.857 198.9357
## alk_cat=(872,1.26e+03] 78 78 78 28 2942.289 201.0198
## alk_cat=(1.26e+03,1.98e+03] 78 78 78 31 2954.372 206.8846
## alk_cat=(1.98e+03,1.39e+04] 78 78 78 48 2454.019 186.7535
## median 0.95LCL 0.95UCL
## alk_cat=[289,872] NA NA NA
## alk_cat=(872,1.26e+03] 3222 2540 NA
## alk_cat=(1.26e+03,1.98e+03] 3574 2769 NA
## alk_cat=(1.98e+03,1.39e+04] 2256 1444 3358
kmast <- survfit(Surv(time, cens) ~ ast_cat, data = daneKNN)
summary(kmast)$table
## records n.max n.start events rmean se(rmean) median
## ast_cat=[26.4,80.6] 79 79 79 18 3721.073 166.1474 NA
## ast_cat=(80.6,115] 78 78 78 24 3299.247 196.2312 3853
## ast_cat=(115,152] 79 79 79 37 2600.008 200.7322 2796
## ast_cat=(152,457] 76 76 76 46 2225.745 197.0123 1690
## 0.95LCL 0.95UCL
## ast_cat=[26.4,80.6] 3762 NA
## ast_cat=(80.6,115] 3282 NA
## ast_cat=(115,152] 2288 3584
## ast_cat=(152,457] 1165 3445
W przypadku zmiennej plat prawdopodobieństwo przeżycia jest większe w grupach, w których pacjenci charakteryzowali się wyższą ilością trombocytów we krwi - funkcje przetrwania dla trzech przedziałów o najwyższych wartościach się pokrywają w pewnym stopniu, natomiast grupa, w której pacjenci mieli najmniejszą ilość płytek krwi, ma najmniejsze prawdopodobieństwo przeżycia. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value < 0.0001). Analogiczna sytuacja występuje w przypadku zmiennej dotyczącej krzepnięcia krwi - im dłuższy czas krzepnięcia krwi, tym prawdopodobieństwo przeżycia jest mniejsze. Test Log-Rank wykazał istotne różnice rozkładów czasu trwania pomiędzy badanymi grupami (p-value < 0.0001).
groupSplots8 <- list()
groupSplots8[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ plat_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"))
groupSplots8[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ pro_cat, data = daneKNN),
pval = TRUE, conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
linetype = "strata",
surv.median.line = "hv",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"))
arrange_ggsurvplots(groupSplots8, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
kmplat <- survfit(Surv(time, cens) ~ plat_cat, data = daneKNN)
summary(kmplat)$table
## records n.max n.start events rmean se(rmean) median
## plat_cat=[62,200] 80 80 80 48 2218.458 180.8581 2105
## plat_cat=(200,256] 77 77 77 27 3151.754 206.5502 3358
## plat_cat=(256,322] 77 77 77 22 3255.308 204.4453 3762
## plat_cat=(322,563] 78 78 78 28 3282.398 186.8769 NA
## 0.95LCL 0.95UCL
## plat_cat=[62,200] 1487 3428
## plat_cat=(200,256] 2503 NA
## plat_cat=(256,322] 3395 NA
## plat_cat=(322,563] 2847 NA
kmpro <- survfit(Surv(time, cens) ~ pro_cat, data = daneKNN)
summary(kmpro)$table
## records n.max n.start events rmean se(rmean) median
## pro_cat=[9,10] 89 89 89 16 3634.561 205.1407 4079
## pro_cat=(10,10.6] 85 85 85 24 3506.679 166.2759 3853
## pro_cat=(10.6,11.1] 61 61 61 31 2801.825 215.3994 3086
## pro_cat=(11.1,17.1] 77 77 77 54 1925.368 195.4251 1191
## 0.95LCL 0.95UCL
## pro_cat=[9,10] 4079 NA
## pro_cat=(10,10.6] 3222 NA
## pro_cat=(10.6,11.1] 2105 NA
## pro_cat=(11.1,17.1] 853 3090
Hazard wzrasta szybciej w grupie mężczyzn niż w grupie kobiet. W przypadku zmiennej stage, największą wartość hazardu mają pacjenci z najwyższym stadium rozwoju choroby.
groupHplots <- list()
groupHplots[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ sex, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_gray(),
palette = c("navajowhite3","plum4"),
fun = "cumhaz")
groupHplots[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ stage, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_minimal(),
palette = c("palegreen3", "skyblue3","salmon3","brown3"),
fun = "cumhaz")
arrange_ggsurvplots(groupHplots, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
Pacjenci, u których zostało stwierdzone powiększenie wątroby (hepato) mają większy hazard niż pacjenci bez tej dolegliwości. Wraz ze wzrostem czasu, różnica hazardów powiększa się. Dla zmiennej skategoryzowanej trig hazard jest najwyższy u pacjentów, którzy posiadają największe wartości trójglicerydów, lecz różnice nie są tak bardzo widoczne.
groupHplots <- list()
groupHplots[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ hepato, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_gray(),
palette = c("palegreen4", "indianred3"),
fun = "cumhaz")
groupHplots[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ trig_cat, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"),
fun = "cumhaz")
arrange_ggsurvplots(groupHplots, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
Hazard jest największy u osób posiadających dużą ilość miedzi w moczu oraz u osób, które posiadają największe stężenie barwnika żółciowego. Pacjenci należący do kategorii z mniejszymi wartościami obydwu zmiennych mają też odpowiednio niższe wartości hazardu.
groupHplots <- list()
groupHplots[[1]] <- ggsurvplot(survfit(Surv(time, cens) ~ copper_cat, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_gray(),
palette = c("wheat3","seagreen3", "slateblue2", "sienna3"),
fun = "cumhaz")
groupHplots[[2]] <- ggsurvplot(survfit(Surv(time, cens) ~ bili_cat, data = daneKNN),
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_minimal(),
palette = c("peachpuff3","lightgreen", "paleturquoise4", "orchid4"),
fun = "cumhaz")
arrange_ggsurvplots(groupHplots, print = TRUE,
ncol = 2, nrow = 1, risk.table.height = 0.2)
Modele regresji dla wybranych fukncji zmiennej losowej T, mają za zadanie pokazać relację między analizowaną funkcją, definiujacą rozkład zmiennej losowej T, a zbiorem zmiennych charakteryzujących jednostki należące do badanej populacji. W analizie przeżycia najczęściej modelowanymi funkcjami są funkcja hazardu, funkcja trwania, a czasami także dystrybuanta. W modelu zaproponowanych przez Coxa [1972], funkcja hazardu bazowego szacowana jest nieparametrycznie dla momentów tk, w których wystąpiło zdarzenie, czyli jest to model semiparametryczny. Model Coxa nazywany jest modelem proporcjonalnych hazardów (PH), co oznacza, że dla dwóch jednostek z wektorami zmiennych x i x*, iloraz ich hazardów nie zależy od czasu t.
Globalny test Walda (p=<2e-16) wskazuje, że przynajmniej 1 zmienna objaśniająca w modelu statystycznie różni sie od 0, czyli jest istotna statystycznie. Zmienne: sex, ascites, spiders, chol, alk.phos, ast, trig oraz platelet nie sa różne od 0, bo p-value > 0,05.
Cox <- coxph(Surv(time,cens)~age+sex+ascites+hepato+spiders+bili+chol+albumin+copper+alk.phos+ast+trig+platelet+protime, data = daneKNN)
summary(Cox)
## Call:
## coxph(formula = Surv(time, cens) ~ age + sex + ascites + hepato +
## spiders + bili + chol + albumin + copper + alk.phos + ast +
## trig + platelet + protime, data = daneKNN)
##
## n= 312, number of events= 125
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 2.706e-02 1.027e+00 1.033e-02 2.619 0.008815 **
## sex0 -1.019e-01 9.031e-01 2.972e-01 -0.343 0.731787
## ascites1 3.428e-01 1.409e+00 3.411e-01 1.005 0.314994
## hepato1 4.667e-01 1.595e+00 2.257e-01 2.068 0.038625 *
## spiders1 1.863e-01 1.205e+00 2.222e-01 0.838 0.401818
## bili 8.917e-02 1.093e+00 2.174e-02 4.102 4.1e-05 ***
## chol 1.604e-04 1.000e+00 3.945e-04 0.407 0.684291
## albumin -1.047e+00 3.511e-01 2.608e-01 -4.013 6.0e-05 ***
## copper 2.949e-03 1.003e+00 1.150e-03 2.564 0.010349 *
## alk.phos -2.074e-05 1.000e+00 3.810e-05 -0.544 0.586184
## ast 3.255e-03 1.003e+00 1.731e-03 1.880 0.060082 .
## trig -1.136e-03 9.989e-01 1.253e-03 -0.906 0.364814
## platelet -3.698e-04 9.996e-01 1.118e-03 -0.331 0.740890
## protime 3.053e-01 1.357e+00 8.253e-02 3.700 0.000216 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0274 0.9733 1.0068 1.0484
## sex0 0.9031 1.1072 0.5044 1.6172
## ascites1 1.4088 0.7098 0.7219 2.7494
## hepato1 1.5947 0.6271 1.0247 2.4818
## spiders1 1.2048 0.8300 0.7794 1.8623
## bili 1.0933 0.9147 1.0477 1.1409
## chol 1.0002 0.9998 0.9994 1.0009
## albumin 0.3511 2.8482 0.2106 0.5854
## copper 1.0030 0.9971 1.0007 1.0052
## alk.phos 1.0000 1.0000 0.9999 1.0001
## ast 1.0033 0.9968 0.9999 1.0067
## trig 0.9989 1.0011 0.9964 1.0013
## platelet 0.9996 1.0004 0.9974 1.0018
## protime 1.3571 0.7369 1.1544 1.5953
##
## Concordance= 0.851 (se = 0.019 )
## Likelihood ratio test= 187.6 on 14 df, p=<2e-16
## Wald test = 210.4 on 14 df, p=<2e-16
## Score (logrank) test = 317.7 on 14 df, p=<2e-16
W wyniku tej procedury otrzymano model z 7 zmiennymi objasniajacymi: age,hepato,bili,albumin,copper,ast,protime. Według kryterium AIC zmienna ast mimo p-value jest równe 0,06 powinna znalezc sie w modelu. Pozostałe zmienne na poziomie istotnosci 0,05 sa statystycznie rozne od 0 (czyli istotne), ponieważ p-value < 0,05.
stepAIC(Cox,direction="backward")
## Start: AIC=1120.33
## Surv(time, cens) ~ age + sex + ascites + hepato + spiders + bili +
## chol + albumin + copper + alk.phos + ast + trig + platelet +
## protime
##
## Df AIC
## - platelet 1 1118.4
## - sex 1 1118.5
## - chol 1 1118.5
## - alk.phos 1 1118.6
## - spiders 1 1119.0
## - trig 1 1119.2
## - ascites 1 1119.3
## <none> 1120.3
## - ast 1 1121.7
## - hepato 1 1122.7
## - copper 1 1124.4
## - age 1 1125.2
## - protime 1 1129.6
## - bili 1 1132.2
## - albumin 1 1134.0
##
## Step: AIC=1118.44
## Surv(time, cens) ~ age + sex + ascites + hepato + spiders + bili +
## chol + albumin + copper + alk.phos + ast + trig + protime
##
## Df AIC
## - chol 1 1116.5
## - sex 1 1116.6
## - alk.phos 1 1116.8
## - spiders 1 1117.2
## - trig 1 1117.5
## - ascites 1 1117.5
## <none> 1118.4
## - ast 1 1120.6
## - hepato 1 1121.2
## - copper 1 1122.4
## - age 1 1123.5
## - protime 1 1127.9
## - bili 1 1130.4
## - albumin 1 1132.3
##
## Step: AIC=1116.53
## Surv(time, cens) ~ age + sex + ascites + hepato + spiders + bili +
## albumin + copper + alk.phos + ast + trig + protime
##
## Df AIC
## - sex 1 1114.7
## - alk.phos 1 1114.9
## - spiders 1 1115.3
## - trig 1 1115.5
## - ascites 1 1115.6
## <none> 1116.5
## - ast 1 1119.0
## - hepato 1 1119.3
## - copper 1 1120.4
## - age 1 1121.5
## - protime 1 1125.9
## - bili 1 1130.2
## - albumin 1 1130.3
##
## Step: AIC=1114.7
## Surv(time, cens) ~ age + ascites + hepato + spiders + bili +
## albumin + copper + alk.phos + ast + trig + protime
##
## Df AIC
## - alk.phos 1 1113.1
## - spiders 1 1113.4
## - trig 1 1113.6
## - ascites 1 1113.9
## <none> 1114.7
## - ast 1 1117.5
## - hepato 1 1117.6
## - copper 1 1120.2
## - age 1 1121.0
## - protime 1 1124.0
## - bili 1 1128.2
## - albumin 1 1128.4
##
## Step: AIC=1113.14
## Surv(time, cens) ~ age + ascites + hepato + spiders + bili +
## albumin + copper + ast + trig + protime
##
## Df AIC
## - trig 1 1112.1
## - spiders 1 1112.2
## - ascites 1 1112.7
## <none> 1113.1
## - hepato 1 1116.0
## - ast 1 1116.1
## - copper 1 1118.2
## - age 1 1120.6
## - protime 1 1122.1
## - albumin 1 1126.4
## - bili 1 1126.5
##
## Step: AIC=1112.11
## Surv(time, cens) ~ age + ascites + hepato + spiders + bili +
## albumin + copper + ast + protime
##
## Df AIC
## - ascites 1 1111.0
## - spiders 1 1111.2
## <none> 1112.1
## - hepato 1 1114.9
## - ast 1 1115.7
## - copper 1 1117.2
## - age 1 1120.9
## - protime 1 1121.9
## - bili 1 1125.2
## - albumin 1 1125.2
##
## Step: AIC=1111.05
## Surv(time, cens) ~ age + hepato + spiders + bili + albumin +
## copper + ast + protime
##
## Df AIC
## - spiders 1 1110.4
## <none> 1111.0
## - hepato 1 1113.6
## - ast 1 1114.2
## - copper 1 1118.5
## - protime 1 1121.7
## - age 1 1121.9
## - bili 1 1126.2
## - albumin 1 1128.5
##
## Step: AIC=1110.39
## Surv(time, cens) ~ age + hepato + bili + albumin + copper + ast +
## protime
##
## Df AIC
## <none> 1110.4
## - ast 1 1113.4
## - hepato 1 1114.0
## - copper 1 1119.1
## - age 1 1120.7
## - protime 1 1122.4
## - bili 1 1127.5
## - albumin 1 1129.1
## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + bili + albumin +
## copper + ast + protime, data = daneKNN)
##
## coef exp(coef) se(coef) z p
## age 0.0326504 1.0331893 0.0092752 3.520 0.000431
## hepato1 0.5080490 1.6620454 0.2174273 2.337 0.019458
## bili 0.0854208 1.0891753 0.0177598 4.810 1.51e-06
## albumin -1.1054157 0.3310732 0.2348614 -4.707 2.52e-06
## copper 0.0032776 1.0032830 0.0009424 3.478 0.000505
## ast 0.0037099 1.0037168 0.0015580 2.381 0.017255
## protime 0.3220291 1.3799250 0.0778588 4.136 3.53e-05
##
## Likelihood ratio test=183.5 on 7 df, p=< 2.2e-16
## n= 312, number of events= 125
Cox1 <- coxph(Surv(time,cens)~age+hepato+bili+albumin+copper+ast+protime, data = daneKNN)
summary(Cox1)
## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + bili + albumin +
## copper + ast + protime, data = daneKNN)
##
## n= 312, number of events= 125
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.0326504 1.0331893 0.0092752 3.520 0.000431 ***
## hepato1 0.5080490 1.6620454 0.2174273 2.337 0.019458 *
## bili 0.0854208 1.0891753 0.0177598 4.810 1.51e-06 ***
## albumin -1.1054157 0.3310732 0.2348614 -4.707 2.52e-06 ***
## copper 0.0032776 1.0032830 0.0009424 3.478 0.000505 ***
## ast 0.0037099 1.0037168 0.0015580 2.381 0.017255 *
## protime 0.3220291 1.3799250 0.0778588 4.136 3.53e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0332 0.9679 1.0146 1.0521
## hepato1 1.6620 0.6017 1.0853 2.5452
## bili 1.0892 0.9181 1.0519 1.1278
## albumin 0.3311 3.0205 0.2089 0.5246
## copper 1.0033 0.9967 1.0014 1.0051
## ast 1.0037 0.9963 1.0007 1.0068
## protime 1.3799 0.7247 1.1846 1.6074
##
## Concordance= 0.845 (se = 0.019 )
## Likelihood ratio test= 183.5 on 7 df, p=<2e-16
## Wald test = 199.3 on 7 df, p=<2e-16
## Score (logrank) test = 276.9 on 7 df, p=<2e-16
Kluczowym założeniem, które musi być spełnione, by model prawidłowo odzwierciedlał wpływ zmiennych objaśniających na rozkład czasu trwania, jest założenie proporcjonalności hazardów. Spełnienie tego założenia powinno być zweryfikowane dla każdej ze zmiennych objaśniających w modelu. Jeżeli p > 0,05 to zostaly spełnione zalożenia proporcjonalności. Dla zmiennych: bili oraz protime te założenia nie zostały spełnione, więc należy usunąć te zmienne z modelu.
Prop <- cox.zph(Cox1, transform = 'km')
print(Prop)
## chisq df p
## age 0.3590 1 0.5491
## hepato 0.8706 1 0.3508
## bili 7.3198 1 0.0068
## albumin 1.1252 1 0.2888
## copper 0.0145 1 0.9043
## ast 1.9955 1 0.1578
## protime 5.2823 1 0.0215
## GLOBAL 16.9008 7 0.0180
par(mfrow=c(4,2))
plot(Prop)
Cox2 <- coxph(Surv(time,cens)~age+hepato+albumin+copper+ast, data = daneKNN)
summary(Cox2)
## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + albumin + copper +
## ast, data = daneKNN)
##
## n= 312, number of events= 125
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.0330674 1.0336202 0.0088058 3.755 0.000173 ***
## hepato1 0.7046633 2.0231653 0.2112452 3.336 0.000851 ***
## albumin -1.2051939 0.2996339 0.2336910 -5.157 2.51e-07 ***
## copper 0.0046068 1.0046175 0.0008553 5.386 7.19e-08 ***
## ast 0.0047647 1.0047761 0.0012846 3.709 0.000208 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0336 0.9675 1.0159 1.0516
## hepato1 2.0232 0.4943 1.3373 3.0609
## albumin 0.2996 3.3374 0.1895 0.4737
## copper 1.0046 0.9954 1.0029 1.0063
## ast 1.0048 0.9952 1.0022 1.0073
##
## Concordance= 0.824 (se = 0.02 )
## Likelihood ratio test= 140.4 on 5 df, p=<2e-16
## Wald test = 162.9 on 5 df, p=<2e-16
## Score (logrank) test = 174.9 on 5 df, p=<2e-16
Metoda skalowanych reszt Schoenfelda Ocenę tego, czy założenie proporcjonalności jest spełnione dokonuje się wzrokowo poprzez sprawdzenie czy punkty na wykresie układają się wzdłuż linii poziomej. Wyraźna zaobserwowana tendencja w przebiegu tej krzywej może sugerować brak proporcjonalności zmiennej.
Ocena na podstawie wykresy wydaje się być problematyczna, ale można zauważyć, że zmienne bill i protime nie spełniają założeń proporcjonalności.
reszty <- resid(Cox1, type="scaledsch")
Time <- as.numeric(rownames(reszty))
zmienne <- names(daneKNN[,c(5,8,11,13,14,16,19)])
par(mfrow = c(3,3))
for (i in 1:7) {
plot(log(Time), reszty[,i], xlab="ln(Czas)", main=zmienne[i],
ylab="Skalowane reszty Schoenfelda", pch=20, cex=0.7)
lines(smooth.spline(log(Time), reszty[,i] ), lwd=3 )
}
Ważnym aspektem ewaluacji modeli Coxa jest badanie jednostek odstających, tj. tych które:
Są to punkty na wykresach oddalone od pozostałych. Jednostki o nietypowej konfiguracji zmiennych objaśniających. Mogą w nieproporcjonalny sposób wpływać na oszacowania parametrów strukturalnych modelu. Do identyifikacji wartości nietypowych wykorzystuje się reszty odchyleń. Za wartości nietypowe można uznać te jednostki, dla których reszty odchyleń przyjmują wartości +- 3 lub więcej.
deviance <- residuals(Cox2,type="deviance")
s <- Cox2$linear.predictors
plot(s,deviance,xlab="Liniowy predyktor",ylab="Reszty odchylen",cex=0.5, pch=20)
abline(h=c(3,-3),lty=3)
daneKNN$deviance <- deviance
c1 <- which(daneKNN$deviance < c(-3))
W przypadku wysokich wartości przyrostu oceny parametru należy rozważyć możliwość usunięcia danej jednostki z próby. Do usunięcia z próby zakwalifikowały się jednostki o numerze 48, 166, 231, 233 i 253. Są one odbiegające i wpływowe.
dfb <- residuals(Cox2,type="dfbeta")
n <- dim(dfb)[1]
obs.nr <- c(1:n)
par(mfrow = c(2,3))
for (j in 1:5) {
plot(obs.nr,dfb[,j],xlab="Numer jednostki",ylab="Przyrost oceny parametru",
main=zmienne[j])
}
a1 <- which(abs(dfb[,1])>(0.004)) #usunac te jednostki
a2 <- which(abs(dfb[,2])>(0.06))
a3 <- which(abs(dfb[,3])>(0.1))
a4 <- which(abs(dfb[,4])>(0.0004))
a5 <- which(abs(dfb[,5])>(0.001))
c <- sort(unique(c(a1,a2,a3,a4,a5,c1)))
daneKNN_1 <- daneKNN[-c,] #zredukowany zbior
Cox3 <- coxph(Surv(time,cens)~age+hepato+albumin+copper+ast, data = daneKNN_1)
summary(Cox3)
## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + albumin + copper +
## ast, data = daneKNN_1)
##
## n= 307, number of events= 124
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.0486020 1.0498024 0.0098801 4.919 8.69e-07 ***
## hepato1 0.6956171 2.0049458 0.2086548 3.334 0.000857 ***
## albumin -1.4274764 0.2399136 0.2545430 -5.608 2.05e-08 ***
## copper 0.0064491 1.0064699 0.0009382 6.874 6.25e-12 ***
## ast 0.0082022 1.0082359 0.0016792 4.885 1.04e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0498 0.9526 1.0297 1.0703
## hepato1 2.0049 0.4988 1.3320 3.0179
## albumin 0.2399 4.1682 0.1457 0.3951
## copper 1.0065 0.9936 1.0046 1.0083
## ast 1.0082 0.9918 1.0049 1.0116
##
## Concordance= 0.836 (se = 0.019 )
## Likelihood ratio test= 180.3 on 5 df, p=<2e-16
## Wald test = 187.2 on 5 df, p=<2e-16
## Score (logrank) test = 217.2 on 5 df, p=<2e-16
Jeżeli zmienna objaśniająca 𝑋 jest ilościowa, to wprowadzenie jej do modelu Coxa jest równoważne przyjęciu założenia, że łączy ją związek liniowy z logarytmem hazardu. Jedna z najcześciej spotykanych metod graficznych jest metoda reszt martyngałowych. Polega na wyznaczeniu reszt martyngałowych dla modelu bez zmiennych objaśniających i zestawieniu ich z badaną zmienną objaśniającą na wykresie korelacyjnym.
Jeżeli funkcja ta ma kształt zbliżony do liniowego to potwierdza to liniowość relacji między zmienną objaśniającą a logarytmem hazardu, w przeciwnym wypadku kształt funkcji może być wskazówką co do wyboru odpowiedniej metody transformacji zmiennej.
age - Potwierdzone założenie o liniowości
zmn1 <- daneKNN_1$age
lab1 <- "Age (years)"
reszty1 <- resid(coxph(Surv(time,cens)~1,data=daneKNN_1),type="martingale")
plot(zmn1, reszty1, xlab=lab1,ylab="Martingale residuals",cex=0.6)
lines(lowess(zmn1, reszty1,delta=1),lwd=2)
albumin - Potwierdzone założenie o liniowości
zmn2 <- daneKNN_1$albumin
lab2 <- "albumin"
reszty2 <- resid(coxph(Surv(time,cens)~1,data=daneKNN_1),type="martingale")
plot(zmn2, reszty2, xlab=lab2,ylab="Reszty martynga?owe",cex=0.6)
lines(lowess(zmn2, reszty2,delta=1),lwd=2)
copper - Brak potwierdzenia założenia o liniowości
zmn3 <- daneKNN_1$copper
lab3 <- "copper"
reszty3 <- resid(coxph(Surv(time,cens)~1,data=daneKNN_1),type="martingale")
plot(zmn3, reszty3, xlab=lab3,ylab="Reszty martynga?owe",cex=0.6)
lines(lowess(zmn3, reszty3,delta=1),lwd=2)
ast - Potwierdzone założenie o liniowości
zmn4 <- daneKNN_1$ast
lab4 <- "ast"
reszty4 <- resid(coxph(Surv(time,cens)~1,data=daneKNN_1),type="martingale")
plot(zmn4, reszty4, xlab=lab4,ylab="Reszty martynga?owe",cex=0.6)
lines(lowess(zmn4, reszty4,delta=1),lwd=2)
W przypadku podejrzenia braku liniowości związku, stosowane są różne metody transformacji zmiennej objaśniającej. Najpopularniejsze z nich, to: dychotomizacja zmiennej objaśniającej (wyznaczenie optymalnego punktu odcięcia)
Dychotomizacja zmiennej ciągłej to zastąpienie jej przez zmienną zero-jedynkową taką, że zmienna ta przyjmuje wartość 0, gdy zmienna ciągła przyjmuje wartości równe lub mniejsze niż punkt odcięcia, oraz wartość 1, gdy zmienna ciągła przyjmuje wartości większe niż punkt odcięcia.
Wartość kryterium informacyjnego Akaike po dychotomizacji okazała się dużo gorsza niż przed tym zabiegiem.
Cox4 <- coxph(Surv(time, cens)~zmn3,data=daneKNN_1)
A1 <- round(AIC(Cox3),2) #akaike z modelu coxa
A1
## [1] 1093.72
cutpoint <- cutp(Cox4)$zmn3 ## optymalny punkt odciecia
c <- (zmn3>=cutpoint$zmn3[1])*1+0 #jezeli zmienna jest >= punkotwi odciecia to przypisuje 1 w przeciwnym wypadku 0
Cox5 <- coxph(Surv(time, cens)~c,data=daneKNN_1)
A2 <- round(AIC(Cox5),2)
A2
## [1] 1201.05
Metoda polega na poszukiwaniu dla zmiennej 𝑋, funkcji 𝑔(𝑥), która najlepiej odzwierciedla prawdziwą relację zmiennej 𝑋 i logarytmu funkcji hazardu.
Wartość kryterium informacyjnego Akaike ponownie jest wyższa w porównaniu do zmiennej copper przed transformacją.
m3 <- mfp(Surv(time,cens)~ fp(zmn3, df = 4, select = 0.05),family=cox, data=daneKNN_1) ## wielomiany ulamkowe
A3 <- round(AIC(m3))
A3
## [1] 1182
int <- quantile(na.omit(zmn3),probs=c(0.05,0.275,.5,.725,.95))
Cox6 <- coxph(Surv(time,cens)~ns(zmn3,knots=int),data=daneKNN_1)
pred <- predict(Cox6,type="terms",se.fit=TRUE, terms=1)
plot(na.omit(zmn3),exp(pred$fit),type="n",xlab=lab3,ylim=c(0,2.5),
ylab="Hazard względny")
lines(smooth.spline(na.omit(zmn3),exp(pred$fit+1.96*pred$se.fit)),lty=2)
lines(smooth.spline(na.omit(zmn3),exp(pred$fit-1.96*pred$se.fit)),lty=2)
lines(smooth.spline(na.omit(zmn3),exp(pred$fit)),lty=1)
abline(h=1,lty=3)
legend('topright',2,c("splajn","95% przedzial ufnosci"), lty=1:2, box.lty=0)
A4 <- round(AIC(Cox6),2)
A4
## [1] 1187.47
Mimo wszystko najlepszym modelem jest model bez transformacji zmiennej copper. Interpretacja parametrów w tym modelu przedstawia się w następujący sposób:
rbind(A1,A2,A3,A4)
## [,1]
## A1 1093.72
## A2 1201.05
## A3 1182.00
## A4 1187.47
Cox3 <- coxph(Surv(time,cens)~age+hepato+albumin+copper+ast, data = daneKNN_1)
summary(Cox3)
## Call:
## coxph(formula = Surv(time, cens) ~ age + hepato + albumin + copper +
## ast, data = daneKNN_1)
##
## n= 307, number of events= 124
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.0486020 1.0498024 0.0098801 4.919 8.69e-07 ***
## hepato1 0.6956171 2.0049458 0.2086548 3.334 0.000857 ***
## albumin -1.4274764 0.2399136 0.2545430 -5.608 2.05e-08 ***
## copper 0.0064491 1.0064699 0.0009382 6.874 6.25e-12 ***
## ast 0.0082022 1.0082359 0.0016792 4.885 1.04e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0498 0.9526 1.0297 1.0703
## hepato1 2.0049 0.4988 1.3320 3.0179
## albumin 0.2399 4.1682 0.1457 0.3951
## copper 1.0065 0.9936 1.0046 1.0083
## ast 1.0082 0.9918 1.0049 1.0116
##
## Concordance= 0.836 (se = 0.019 )
## Likelihood ratio test= 180.3 on 5 df, p=<2e-16
## Wald test = 187.2 on 5 df, p=<2e-16
## Score (logrank) test = 217.2 on 5 df, p=<2e-16
Poniżej ukazano miary dopasowania R2, (obliczone na podstawie statystyki cząstkowej ilorazu wiarygodności w modelu Coxa) Najlepiej dopasowanym modelem (o największej wartości tej miary = 0,77) jest model Cox3 otrzymany metodą krokową, po odrzuceniu zmiennych niespełniajacych założeń proporcjonalnosci oraz jednostek odstających i wpływowych. Co ciekawe, taką samą wartość miary R2 (0,77) obliczono dla modelu Cox1 otrzymanego metodą krokową przed sprawdzeniem zalożeń proporcjonalnosci oraz przed odrzuceniem jednostek odstających i wpływowych
r2_1 <- coxr2(Cox1) #model otrzymany metoda krokowa przed sprawdzeniem zalozen proporcjonalnosci
r2_2 <- coxr2(Cox2) #model po odrzuceniu zmiennych niespelniajacych zalozen proporcjonalnosci
r2_3 <- coxr2(Cox3) #model po odrzuceniu jednostek odstajacych i wplywowych (uznany za najlepszy)
round(rbind(r2_1$rsq, r2_2$rsq, r2_3$rsq),2)
## rsq
## [1,] 0.77
## [2,] 0.67
## [3,] 0.77